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Abstract 

We describe a stochastic technique which allows one to compute numeri- 
cally the coefficients of the weak coupling perturbative expansion of any 
observable in Lattice Gauge Theory. The idea is to insert the exponen- 
tial representation of the link variables U^(x) — > exp{A^(x) / 'vT?} into the 

Langevin algorithm and the observables and to perform the expansion in 
_i 

(3 2. The Langevin algorithm is converted into an infinite hierarchy of 
maps which can be exactly truncated at any order. We give the result 
for the simple plaquette of SU(3) up to fourth loop order (/3 -4 ) which 
extends by one loop the previously known series. 
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1 Introduction 



Weak coupling perturbation theory has been developed since the early times of Lat- 
tice Gauge Theory While the lattice formulation provides the only method which 
allows for systematic non-perturbative calculations from first principles in strongly 
interacting asymptotically free theories, its perturbative expansion is also important. 
In particular weak coupling perturbation theory on the lattice is useful in order to 
accelerate the approach to the continuum, like in Symanzik's improving program, or 
in the calculation of renormalization constants. For observables such as the gluon 
condensate or the topological susceptibility the weak coupling perturbative expan- 
sion on the lattice gives an important contribution to the Monte Carlo data (see 
Ref. |2], [| [|). The expansion coefficients of the lattice correction to the gluon con- 
densate and the topological susceptibility have been computed to three loop order || 
by means of a rather sophisticated diagrammatic technique. 

We recently proposed || an alternative technique to compute numerically the 
coefficients of the expansion of any lattice observable based on Parisi-Wu J!]] stochas- 
tic quantization implemented on the lattice ||. The idea is simple: we insert the 
exponential map U^(x) = exp{A /i (x) / V/?} everywhere in the elementary Langevin 
evolution step. The stochastic evolution of the Lie-algebraic field depends para- 
metrically on (3. We expand in powers of an d obtain the evolution in the 

Langevin time order by order. By expanding the observable under consideration, 
the plaquette variable in the present study, it turns out that the coefficients of the 
perturbative expansion are given by expectation values of composite operators and 
are estimated as an average on the Langevin history. 

A fundamental aspect of the approach is concerned with the problem of gauge 
fixing. We have found that the Langevin evolution with no gauge fixing is affected 
by a random fluctuation increasing in time and that this problem is taken care of 
by adopting a stochastic gauge fixing idea proposed by Zwanziger |J and already 
implemented on the lattice |10| . 



In this paper we present an application of the method by computing the four loop 
coefficient of the plaquette expectation value which is related to the expectation value 
of the condensate ^(F 2 (x)}. The first three coefficients agree with the known values. 

We present the algorithm in some detail in section two. In section three we recall 
the stochastic gauge fixing method which is needed to keep statistical fluctuations 
finite. In section four we report the results obtained for the simple plaquette up 
to fourth loop order. In the final section we discuss the possibility of extracting 
the condensate from the plaquette Monte Carlo data and we comment on a possible 
scenario of high order perturbation theory. 
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2 The algorithm 



We consider the standard pure SU (N) gauge lattice action 

^] = -^£Tr(f/ P + 4) (1) 

where the sum extends to all plaquettes P in a hypercubic lattice. We shall denote by 
Ufj,(x) the link variable at the site x in the direction \i. We borrow the implementation 
of the Langevin dynamics from Ref. ||. A single Langevin step is given by a sweep of 
the lattice where each link variable is updated according to the rule 

U,(x) -> %{x) = U,(x) (2) 

where F^(x) is given by 

W = 7§ E VP ~ U P )\traceless + ^H^x) (3) 

Here e is the Langevin time step; the sum over P means that gets contributions 
from all oriented plaquettes which include the link fi at x\ the suffix traceless means 
that a subtraction of l/iVTr(.) times the identity matrix is understood. Finally 
H^ix) is extracted from a standard (antihermitian, traceless) Gaussian ensemble of 
A-dimensional matrices. Now we represent by the exponential map on SU(N): 
U^(x) = exp{A ll (x) I • The Langevin evolution, when translated in terms of A, is 
going to depend parametrically on (3. We then insert a formal power series expansion 

^(^E^) (4) 

fc>0 

into Eq. (^) and we try to match both sides order by order. To make this possible 
it is immediately evident that the time step must be scaled by putting e = r//3, 
which is also expected since going to large (3 requires smaller and smaller e to avoid 
systematic errors. The Langevin algorithm has now been transformed into an infinite 
system of coupled non-linear stochastic finite-difference equations. The system can 
be truncated exactly at the desired order since each field A^jy is independent from 
higher order fields A^ 1 ' with h > k; in particular the field A^ transforms by itself 
and represents the Abelian limit (a collection of (A^ 2 — 1) U(l) gauge fields). The 
explicit form of the Langevin algorithm at any order can be obtained by expanding 
the products of link variables by Baker-Haussdorff-Campbell's formula. Since our 
goal is to reach and, possibly, go beyond the order /3~ 4 , we faced the problem of 
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automatically generating these expansions. This was achieved using Mathematical, 
which generates the algebraic expansion, optimizes the code trying to reduce the 
number of operations to a minimum and finally translates the output to (parallel) 
fortran or apese. Presently we are running the code designed to reach 0(/3~~ 5 ). 

The explicit form of the expansion F^(x) rapidly gets very cumbersome; to get 
a rough idea, the number of lines of text in the Mathematica output starts at 47 at 
order (3~ 2 and it rises to 3544 at order f3~ A . Hence we report the stochastic equations 
only for the first few orders (already specialized to SU(3)) 

Af' = Af-F® 

Af' = Af-F^-l[F^,Af) 

Af = AV-FP-1[F?\AV]-$[FP\A®] 

Af = Af - if) - Af] - |[if \ Af) - i[if \ Af) 

+ [if), Af}] + 1[F(°), [F«, Af}] + l[if>, [if >, Af}} 

~ Y 2 [Af, [Af, if)]] - 1[4°), [Af, if)]] - ±[Af, [Af, if)]] 

- ^W 0) ) [4 0) .[^ 0) .4 0) ]]] ) ( 5 ) 

where we denote by if )(x) the expansion coefficients of the unitary drift F^(x) = 
J2k>o P~^ k+1 ^ 2 Fu k \ x ) w hich are explicitly given to the first two orders by 

= I E Af{x) + ^H,{x) 
Ff\x) = - £ Af(x) + ^- £ [4 0) (^),4 0) (^)] (6) 

v<p 

Given the stochastic evolution order by order, we consider the expansion of any 
observable O ~ f3~ k ^ 2 O k and we get 



(O) = £/r fe / 2 <(9 fc > 



i T 
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where O k is the value of the random variable, averaged over the whole lattice, after t 
steps. The existence of the limit in terms of some "perturbative probability distribu- 
tion" is still to be proven and it will be assumed as a working hypothesis here. This 
fact deserves further consideration in view of the fact that at first sight the behavior 
of the stochastic process is far from trivial. 



3 Stochastic Gauge Fixing 

The first run up to level (3~ 2 was done without any gauge fixing. The evidence 
(see 0) was rather strong in favour of the existence of the time average with an 
excellent agreement with the known values; however the statistical fluctuation of 
the signal at order 1 and 2 tends to increase in time signaling the fact that the 
asymptotic distribution, if it exists at all, has infinite second moments. The danger 
of such wild fluctuations in a numerical calculation is rather obvious - as soon as the 
signal/fluctuation ratio is smaller than the machine accuracy the algorithm is dead. 
While we are not able at present to suggest a precise picture of these fluctuations, 
it is believed that fixing the gauge is also going to fix the problem. Actually these 
large fluctuations can be thought of as gauge-non-invariant contributions to the time 
trajectory (with vanishing expectation) due to the fact that the full gauge invariance is 
only implemented by averaging on all possible noise-histories; for a given noise history 
stochastic time evolution and gauge transformations do not commute and this fact was 



considered responsible for this kind of fluctuations in Langevin dynamics (see ||11|| ). 
The immediate source of divergence can be identified in the fields A^\ which fluctuate 
like free Brownian variables in the absence of gauge-fixing. Higher order fields 
contain powers of and hence show stronger fluctuations. Fixing the gauge keeps 
the fields A^ to bounded deviations and hence cures the problem. The kind of 
gauge fixing appropriate to Langevin dynamics was introduced by Zwanziger ([[J) 
and implemented as a lattice algorithm in [|10J. The method consists in performing 
a gauge transformation after each sweep, in such a way that the field is attracted to 
the manifold defined by Landau gauge. The gauge transformation is given by 

W L : U^x) -> e w{x) U^x)e~ w(x+e ^ 

w(x) = a ^-ni u n( x ) ~ Ul{x)]\\ traC ei ess 

A^U u (x) = U v (x) - U v {x - e„) (7) 

An iteration of Wl by itself would bring the gauge potential to the Landau gauge. 
By interleaving it to each Langevin step one obtains a sort of soft gauge fixing, 
Wl providing an additional drift which however does not modify the asymptotic 
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n 


1 


2 


3 


4 


c n (r = 0.02) 


2.019 


1.210 


3.001 


9.58 


Ac n 


0.001 


0.002 


0.014 


0.09 


c n (r -> 0) 


2.000(2) 


1.218(7) 






Exact 


2 


1.212(7) 


3.12 





Table 1: Stochastic estimate of perturbative coefficients for the simple plaquette, 8 4 
lattice. 

probability distribution. We have now to expand the gauge fixing step Wl in 
this can be achieved with the same technique already developed for the unconstrained 
Lengevin algorithm. The parameter a is chosen in order to minimize systematic errors 
— in practice it is set equal to r. 

4 Results 

All the data will refer to the SU (3) simple plaquette variable 

W=l- i(Tr(C/ P )) ~ ]T/r"/ 2 <a> = E c «/?" n • (8) 

n 

The expansion coefficients are obtained from the time average taken on a trajectory 
of 2100 Langevin steps of which the first 500 were discarded (see Fig. 1). The odd- 
degree operators 02n+i, which must have vanishing expectation, are used to monitor 
the process (see Fig. 2). 

The systematic error due to the finite evolution step can be corrected by running 
at different values of r and extrapolating at r = (see Tab. 1). We check in this way 
that the method gives very accurate results when compared to the analytic results, but 
already for c 2 we need more statistics to get an improvement from the extrapolation 
in r (see Fig. 3). 

Let us now discuss to what extent our calculation of the perturbative coefficients 
is different from those previously published. First of all, of course, there is an inherent 
statistical error, due to the stochastic nature of the calculation, which can be improved 
by accumulating more data. 

Other sources of systematic errors, which are not seen with the present accuracy, 
may turn out to be relevant for other observables. First of all there are finite size 
effects: since we have to subtract the perturbative tail from Monte Carlo data taken 
on the same finite lattice, it could turn out that the right thing to do is not to correct 
for finite size at all. For a local quantity like the gluon condensate the finite size effect 
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is going to be negligible anyway. Another difference from the standard calculation is 
more fundamental: we are not attempting to extract the leading term in the lattice 
spacing. Our result includes all irrelevant higher order terms in a, which are also 
present in Monte Carlo data. 

5 Discussion 

Let us now recall from Ref. H that the single plaquette expectation Wmc obtained 
by Monte Carlo simulations on the lattice is related to the physical condensate by 

W MC = W peTt + - (F 2 ) a 4 Z(f3) , W pcit = £ c n (3~ n , (9) 

77 n>0 

where W pert is the lattice perturbative correction discussed before. To two loop order 
the condensate in SU(3) is given by the asymptotic freedom expression 

G?(J3) = -(F 2 ) a 4 Z(J3) ~ C as (/3/66 ) 26l/6 °e-^ 36 ° (1 - ^ + 0{l/(3 2 )) (10) 

where the constant C as is given by the physical quantity C as = ^(F 2 )/A 4 , with A 
the fundamental lattice gauge theory scale and where 

&o = ll/(47r) 2 , h = 102/(4tt) 4 

are the first two perturbative coefficient of the /3-function. 

The main question is then how to extract the physical condensate from the Monte 
Carlo data. First of all we perform the same analysis done in Ref. we try to 
estimate the gluon condensate by comparing G^ifl) with the various approximations 
given by 

n 

G { ?\(3) = Wmc - W$ t , W$ t = £ c fc /T fc . (11) 

k=l 

In Fig. 4 we report (/3) for n — 1, ... ,4 by using the values of c n in Tab. 1 and 
the Monte Carlo results on the single plaquette reported in Ref. M. By increasing 
the perturbative order n one sees that G^\l3) has a logarithmic slope in (3 which 
decreases monotonically and seems to approach the asymptotic freedom slope in (|I"0"D 
given by the dotted line. The series seems still too short for n — 4. Following the same 
procedure as in Ref. U, we introduce the next two coefficients c 5 and c 6 which we 
obtain by fitting G^\(3) with the asymptotic freedom expression G% s ([3) in Eq. ([10|). 
Together with the two coefficients c 5 , c 6 one fits also the value C as of the condensate. 
The resulting points of G^\/3) are the lower points in Fig. 4 which are nicely fitted 
by the dotted line, the asymptotic freedom expression G^ifi). 
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From the fit one obtains: c 5 = -84, c 6 = 1193, and C as = 3.6 x 10 s . This 
value of C as is close to the estimate obtained by performing the same analysis at one 
lower loop n — 3. Notice that, while the computed coefficients are all positive, the 
first unknown coefficient C5 turns out to be negative. This is simply due to the fact 
that the functions \b.G\ \(/3) have a slight curvature. Indeed, if one applies the same 
procedure at a lower loop, the first fitted coefficient turns out to be negative although 
the known value is positive. 

To analyze more closely the possible approach to asymptotic freedom, we plot in 
Fig. 5 the effective slopes of In G^^ff) computed as follows. We perform a linear fit to 
In G^ifi) in various intervals (fa, . . . , fa) included in the interval 6 < (3 < 7 where we 
expect to observe the asymptotic scaling. The result is shown in Fig. 5, by plotting 
R = 1/ slope vs. 1/y/n, the error bars being a measure of the different estimate of 
the slope depending on the interval chosen; the asymptotic freedom result R as = 3&o 
is plotted as a horizontal line for comparison. It appears that the values of have 
a smooth behaviour as function of l/y/n, which would however extrapolate to R as 
for a finite n. This fact is typical of a divergent series. Indeed it has been observed 



12| , |13l that this should be the case due to the presence of renormalons. The leading 



renormalon gives the following large order behaviour of the expansion coefficients 

c n ~ An B n\ (3b ) n . (12) 



As discussed in Ref. [jTB), this behaviour would imply that the intrinsic error in the 
estimate based on the divergent series would be of the same order of magnitude as the 
asymptotic freedom result Actually for n < 4 we find that the coefficients 

increase even faster with n than the one in Eq. (0). This fact would make even 
more ambiguous the attempt to extract in this manner the gluon condensate from 
the Monte Carlo data. Thus one should find some new criterion to identify the 
perturbative subtraction. The stochastic calculation of c 5 should be available soon 
and we hope it will help in clarifying the situation. 

We conclude by observing that our method can be easily adapted to evaluate high 
order perturbative expansions for any other gauge field observable, in particular we 
may correct for perturbative contributions to Wilson loops and to the topological 
susceptibility; the application to other models, such as chiral models on the lattice, 
should also be straightforward. 
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Figure 1: Time history of 2n , r = 0.02 . 
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Figure 2: Time history of 2n +i-T = 0.02 . 
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Figure 4: at various loop orders. Dashed line is the asymptotic freedom result 

(Eq. (0)) with the constant C as obtained by fitting c 5 and c 6 . The lowest points 
correspond to G% with the fitted coefficients C5 and Cq. The error bars refer to the 
Monte Carlo statistical error only. 
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